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ABSTRACT 

Recently, we have presented the first, truly large-scale radiative transfer simulations of in- 
homogeneous cosmic reionization which resolve all the possible halo sources down to the 
dwarf galaxy mass range, M ^ 10 9 M Q , in a comoving volume (100 /i _1 Mpc) 3 . This is large 
enough to sample the global mean history, geometry and statistical properties of reionization 
fairly and accurately for the first time. Here we present new simulations which extend the 
source halo mass range downward to 10 8 Mq, to capture the full range of halo masses thought 
to be primarily responsible for reionization by their star formation following atomic hydrogen 
radiative cooling and gravitational collapse. Haloes below about 10 9 M Q , however, are sub- 
ject to Jeans-mass filtering in the ionized regions, which suppresses their baryonic content and 
their ability to release ionizing radiation. By including these smaller-mass haloes but account- 
ing for their suppression, too, we find that reionization is "self-regulating," as follows. As 
the mean ionized fraction rises, so does the fraction of the volume within which suppression 
occurs. Hence, the degree of suppression is related to the mean ionized fraction. Since low- 
mass haloes with high efficiency (i.e. high emissivity) achieve a given mean ionized fraction 
earlier than do those with low efficiency, Jeans-mass filtering compensates for the difference 
in the emissivity of the suppressible haloes in these two cases. As a result, in the presence of 
lower-mass source haloes, reionization begins earlier, but the later stages of reionization and 
the time of overlap are dictated by the efficiency of the higher-mass haloes, independent of the 
efficiency of the suppressible, lower-mass haloes. Hence, while the lower-mass haloes do not 
alter the overlap redshift, z ov , they serve to boost the electron-scattering optical depth of the 
universe, r es . This may explain why observations of quasar absorption spectra at high redshift 
find that reionization ended late (z ov < 7), while CMB polarization measurements report a 
large enough r cs that reionization must have begun much earlier (z > 11). We present results 
for the ACDM universe with cosmological parameters from both 1-year and 3-year data re- 
leases of WMAP. Reionization histories consistent with current constraints on z ov and r cs are 
shown to be achievable with standard stellar sources in haloes above 10 8 Af Q . Neither miniha- 
los nor exotic sources are required, and the phenomenon of "double reionization" previously 
suggested does not occur. 

Key words: cosmology: theory — radiative transfer — intergalactic medium — large-scale 
structure of universe — galaxies: formation — radio lines: galaxies 



1 INTRODUCTION 

The inhomogeneous reionization of the intergalactic medium 
(IGM) at high redshift proceeded by the propagation of ioniza- 
tion fronts (I-fronts) outward from the early galaxies that formed 
sources of ionizing radiation, like stars and mini-quasars. This pro- 
cess continued until the H II regions bo unded by these I-fronts 
grew to overlap, fshapiro & Girouxl J 1987h showed that the I-fronts 
during reionization were generally weak, R-type, which move su- 
personically relative to both the neutral gas ahead of and the 
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ionized gas behind them, outracing the hydrodynamic response 
of the IGM to the radiation. This property makes it possible to 
simulate inhomogeneous reionization in a cosmological universe 
like ACDM by performing a radiative transfer calculation on a 
pre-computed, 3-D, cosmological density field generated by ei- 
ther a pure N-body or a gas and N-body dynamics simulation of 
large-scale structure formation. In this approximation - the "static 
limit" - the dark matter and baryonic gas evolve just as they 
would have without reionization, with no back-reaction from the 
gas pressure forces associated with photoheating in an inhomo- 
geneous density field. These structure formation simulations also 
yield the mass and location of the galactic halos which are the 
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sites of star formation and, hence, of ionizing photon release. At- 
temp ts to simulate reionization in a CDM universe this way in- 
cludelCiardi et al.l(l200oh:lNakamoto et al.ld200lh;lRazoumov et alJ 
2003); ICiardi et al.l j2003bh; ISokasian et al.1 ( I2003h : llliev etail 
2006b); Melle maet al.l j2006b| ) . Of theseTthe simul ations de- 



scribed in llliev et alJ d2006b) and Mellema et al. ( 2006b) represent 
the first truly large-scale radiative transfer simulations of reioniza- 
tion, in a comoving volume of (100 h -1 Mpc) 3 , which resolve all 
the individual halo sources down to the dwarf galaxy range, above 
about a billion solar masses. This represents a more than two orders 
of magnitude volume increase over previous work, as is required 
in order to make statistically meaningful predictions of observable 
consequences and features of reionization. To accomplish this, a 
new and efficient radiative transfer method, called C 2 -Ray (Conser - 
vative, Causal Ray-Tracing) was developed (Mell ema et al.l2006al) . 
It was then coupled to the results of very large N-body simulations 
of the ACDM universe involving 1624 3 = 4.28 billion particles 
on a cubic lattice of 3248 3 cells by the parallel PMFAST method 
dMerz et alj|2005l) . 

The back-reaction of the IGM and of the halos which col- 
lapse out of it to make sources cannot be entirely neglected, how- 
ever. These effects bec ome important at sm all scales, in particu- 
lar. As first described in lShapiro et alJ i 19941 henceforth, SGB94), 
the heating of the IGM during reionization introduces pressure 
forces which oppose the growth of linear perturbations in the 
baryonic component by gravitational instability - sometimes re- 
ferred to as "Jeans-mass filtering" - resulting in a negative feed- 
back on the rate of collapse of new baryons out of the IGM into 
dark matter halos. Since small-mass halos form first and merge 
to make larger-mass halos later during hierarchical structure for- 
mation in the ACDM universe, this feedback effect (which pre- 
vents the smaller-mass halo sources from forming inside H II re- 
gions) acts as a "self-regulator" on reionization, as follows. The 
more sources that form, the higher the mean ionized fraction of 
the universe, but the higher the fraction of the IGM which is ion- 
ized, the more of these small-mass halos are suppressed. This slows 
the rapid exponential rise of the source population and, with it, the 
rise of the ionized fraction. Eventually, when halos massive enough 
to overcome this feedback effect and capture their fair share of 
baryons even inside H II regions become common enough, they 
are able to finish reionization without being inhibited by radia- 
tive feedback. This phenomenon of Jeans-mass filtering was subse- 
quently revisited and confirmed by comparison wit h simulations by 
iGnedin & Huil ( 1998) and Gnedinl §000), see also Bi & D avidsenl 
< 19971) . Estimates of the minimum halo mass scale which is suffi- 
cient to overcome the effects of uniform photoionization and ac- 
crete or retain enough baryons to make star formation possible 
have also been refine d over the years, by 1-D and 3-D gasdynam- 
ical simulations (e.g. Efsta thiou Il992l; iThoul & Weinberg! 1 19961 ; 
iNavarro & Steinmetz|[l997l ; iDiikstra et alj|2004l) . Roughly speak- 
ing, for halos during the epoch of reionization (EOR) at redshifts 
z > 6, the minimum halo mass required to avoid suppression is 
about a billion solar masses. 

This mass scale corresponds roughly to the minimum mass 
of the source halos in cluded in t he larg e-scal e reionization 
simulat ions described in llliev et al] d2006bh and iMellema et al.l 
(2006bj) for comovin g simula ti on vol umes of (100 h _1 Mpc) 3 , 
those of ICiardi et al.l d2003bQ 2006), for sim ulation volumes 



(20h _1 Mpc) 3 and those of IZahn et all d2006h for a volume of 
(65 h _1 Mpc) 3 . For these simulations the Jeans-mass filtering is 
not taken into account, since the mass range of the source halos 



they resolved is above the range which is suppressed by reioniza- 
tion. 

The minimum mass of the halos responsible for reionization 
may have been smaller than this, however. The star-forming abil- 
ities of halos are generally thought to be different for halos with 
virial temperatures T v i r above and below 10 4 K, respectively. Halos 
with T vir < 10 4 K (roughly 10 4 < M/M Q < 10 8 ) called "mini- 
haloes" only form stars if they can form a trace of H2 molecules 
sufficient to cool the gas well below T v i r by collisional excita- 
tion of rotational-vibrational lines. Since minihaloes are the ear- 
liest halos to form in the CDM universe, they are expected to 
be the first sites of star formation. Minihaloes are highly vulner- 
able to radiative feedback from these very stars, however. Sim- 
ulations currently suggest that these first stars were massive and 
hot and, henc e, were strong emitters of ionizing an d dissociating 
UV radiation dAbel et al.l 2002; Brom m et al. 1 120021) . A single star 
was then sufficient to expel t he gas from the minihalo that made 
it, by photoionizing the ga s dWhalen et al.l |2004| ; iKitavama et al.l 
|2004| ; lAlvarez et aUl2006af) . The UV background from these stars 
in the Lyman-Werner bands is likely to have dissociated the H2 
molecules in minihaloes which had not yet formed stars, long be- 
fore t hey released enough ionizing photons to reionize the uni- 
verse (Hai manetalE oOO). This suggests that minihaloes were not 
the primary source of reionization, although this conclusion re- 
mains highly uncertain. Those minihaloes which were "sterilized" 
in this way against their own star formation would have trapped 
the intergalactic I-fronts which encountered them during reioniza- 
tion, transforming these I-fronts from R-typ e to D-type, and ex- 
pelling their gas in a photoevaporative wind dShapiro et alj|2004t 
llliev et al.l20 05b). This process might have affected the progress of 
reioniz at ion by consuming additional ionizing photons dlliev et all 
2005bB ICiardi et al.l l2006). The effect of minihaloes as sinks of 
ionizing photons tends to be degenerate, in fact, with the unknown 
efficiency for release of ionizing photons by higher-mass sources, 
as a result of the tendency for minihaloes to cluster around them 
dlliev et all2005ah . 

By contrast with the minihaloes, halos with T v i r > 10 4 K 
(M £ 1O 8 M0) were less vulnerable to the suppression of their star 
formation. These halos were able to cool radiatively and collapse 
by collisional excitation of the atomic hydrogen Lya line, leading 
to star formation even in the presence of the rising UV dissociation 
background from other stars. As described above, however, these 
halos were nevertheless subject to Jeans-mass filtering inside H II 
regions. 

Halos in the intermediate mass range (10 s < M/Mq < 10 9 ), 
therefore, can be important sources of additional ionizing photons 
not explicitly accounted for in our previous large-scale simula- 
tions of reionization. These sources would have been suppressed 
inside H II regions, however, so their overall impact remains to 
be determined. Semi-analytical studies suggest that the inclusion 
of halos subject to Jeans-mass filtering has a significant effect on 
the mean reionization h istory dShapiro et alJI 19941; Ishapirdl 19951; 
I Chiu & Ostrikerl |20QO|; iHaiman & Holderl 120031; Wvithe & Loebl 
120031 ; lOnken & Miralda-Escudel |2004| ; iFurlanetto & Loebl I20051T 
This may help explain why the CMB polarization experi- 
ments indicate that re ionization was well advanced by z ^ 11 
( Spergel & et al. 2006), while quasar absorption spectra suggest 
that reionization ended later, at 2 ~ 6.5 (e.g. iFan et al' 2002; 
IWhite et alj|2003l) , implying that the EOR was extended in time. 

Some smaller-scale reionization simulations have been per- 
formed that couple gas and N-body dynamics to radiative transfer, 
which, in principle, takes account of the feedback effect of Jeans- 
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mass filtering jRicotti et al.ll2002h . However, these simulation vol- 
umes are too small to describe the global reionization history or its 
geometry and statistical properties. As a result, they failed to antic- 
ipate the extended nature of reionization which self-regulation by 
Jeans-mass filtering makes possiblf]That will be the subject of the 
current paper. 

In this paper we assume a flat (flk = 0) ACDM 
cosmology with parameters (f2 m , Ha, flt, h , as, n ) = 
(0.27,0.73,0.044,0.7,0.9,1) fepergel & et all |2003h . here- 
after WMAP1, where f2 m , Qa, and fib are the total matter, 
vacuum, and baryonic densities in units of the critical den- 
sity, as is the rms density fluctuations extrapolated to the 
present on the scale of 8ft -1 Mpc according to the linear per- 
turbation theory, and n is the index of the primordial power 
spectrum of density fluctuations. The ne w, 3 -year WMAP dat a 
yielded somewhat different parameters dSpergel & etail I2006T) . 
(£l m ,tl A ,£l b ,h,a s ,n) = (0.24,0.76,0.042,0.73,0.74,0.95), 
hereafter WMAP3. 

The differences between these WMAP1 and WMAP3 param- 
eters can have a significa nt impact on the progress of reioniza- 
tion d Alvarez et al J 1200 6b'). Structure formation is delayed in the 
WMAP3 universe relative to the WMAP1 universe, especially on 
the small scales relevant to the formation of reionization source 
halos at high redshift, so the epoch of reionization is shifted to 
lower redshifts. In particular, if source halos of a given mass are as- 
sumed to have released ionizing photons with the same efficiency 
in either case, then reionization for WMAP3 is predicted to have 
occurred at (1 + z)-values which are roughly 1.4 times smaller 
than for WMAP1. As such, the predicted electron-scattering opti- 
cal depth of the IGM accumulated since the beginning of the EOR 
would have been smaller for WMAP3 than for WMAP1 by a factor 
of 1.4 3/2 ~ 1.7, just as the observations of large-angle fluctua- 
tions in the CMB polarization require. This means that the ioniz- 
ing efficiency per collapsed baryon required to make reionization 
early enough to explain the value of r cs reported for WMAP1 and 
WMAP3 are nearly the same, despite the fact that r os is smaller for 
WMAP3 than for WMAP1. 

In Section 2, we describe our simulations, the method, input 
parameters and cases. We motivate these simulations and antici- 
pate some of the trends by a simple analytical toy model for the 
mean ionization history in Section 3. Simulation results are pre- 
sented in Section 4, with conclusions in Section 5. Appendix A 
describes how we can use the simulation results to improve upon 
the toy model introduced in Section 3. 



2 THE SIMULATIONS 

Our ba sic methodology was described in d e tail in llliev et all 
(2006b) (hereafter Paper I) and lMellema et alj d2006bt) (hereafter 
Paper II). Hence, here we will outline the main simulation param- 
eters and concentrate on the simulations which have not been pre- 
sented before and the new features we introduce. 

We use very high resolution N-body simulations to derive halo 
catalogues, which include the detailed halo properties, as well as 
their positions and velocities at a number of time-slices (between 50 
and 100 per simulation) and the corresponding gas density fields. 

1 Recently, small-scale reionization simulations like these were also used 
as the basis for a phenomenological "subgrid model" in a large-scale reion- 
ization simulation which was too c oarse-grained to re solve the small-scale 
structure and source halos directly iKohler et al.l200j) . 



All the halos found in the simulation volume at a given redshift 
are assumed to be sources of ionizing radiation (unless they are 
suppressed by Jeans mass filtering, see below). This results typi- 
cally in tens to hundreds of thousands of sources, depending on the 
case. We use a simple recipe to assign a photon emissivity to each 
source, by assuming a constant mass-to-light ratio. We follow the 
time-dependent propagation of the ionization fronts produced by all 
sources in the simulation volume using our detailed radiative trans- 
fer an d non-equilibrium chemistry code C 2 -Ray dMellema et al.l 
l2006ah . which h as been extensively t ested against available ana- 
lytical solutions dMellema et ail2 006a) and a number of other cos- 
mological radiative transfer codes ( llliev et al . 2006i). 

The underlying N-body simulations have a spatial grid of 
3248 3 cells and follow the evolution of 162 4 3 = 4.3 billion parti- 
cles using the particle-mesh code PMFAST dMerz et alj2005h . The 
number of resolved halos for the simulations with WMAP1 cos- 
mology parameters are ~ 4 x 10 5 (~ 8 x 10 5 ) halos at z — 8 
(z = 6)inthe(100h~ 1 Mpc) 3 volume, and ~ 7xl0 5 (~ 9xl0 5 ) 
halos at z = 8 (z = 6) in the (35 h" 1 Mpc) 3 volume. The corre- 
sponding numbers for the WMAP3 simulations are ~ 7.5 x 10 4 
(~3x 10 5 ) halos at z = 8 (z = 6) in the (100 h _1 Mpc) 3 vol- 
ume and ~ 2 x 10 5 (~3x 10 5 ) halos at z = 8 (z = 6) in the 
(35 h" 1 Mpc) 3 volume. 

Simulating the transfer of ionizing radiation with the same 
grid resolution as the underlying N-body is still not feasible on cur- 
rent computers. We therefore re-grid the data to lower resolution, 
with either 203 3 or 406 3 cells, for the radiative transfer simulations. 
We combine sources which fall into the same coarse cell, which re- 
duces slightly the number of sources to be considered compared to 
the total number of halos. 



2.1 Source suppression by Jeans-mass filtering 

The process of photoionization also heats the gas to temperatures 
above 10 4 K. The exact value of the temperature reached varies 
and generally de pends on the local level of the ionizing flux and 
its spectrum (see Shapiro et al.ll2004l for detailed numerical calcu- 
lations). Typical values are Tigm = 10, 000 - 20, 000 K, but it 
could be as high as ~ 40, 000 K for hot (Pop. Ill) black-body spec- 
trum. However, as was mentioned above, the hydrogen line cooling 
is highly efficient for T > 10 4 K, particularly at high redshifts, 
where the gas is denser on average, which would typically bring 
its temperature down to Tigm ~ 10 4 K, and possibly somewhat 
below that due to the adiabatic cooling from the expansion of the 
universe. 

This increase of the IGM temperature caused by its photo- 
heating results in a corresponding increase in the Jeans mass. The 
adiabatic IGM temperature at z ,$ 130, after Compton scattering 
ceases to couple Tigm to T"cmb as it does at earlier times, is well- 
approximated by 

Tigm.o = 26mK(l + z) 2 (1) 

dCouchmanlll985l : lGlover & Branoj|2003l) . In linear theory, the in- 
stantaneous cosmological Jeans mass of the neutral IGM in the ab- 
sence of heating is then given by 

(e.g. lShapiro et alJl994lDiev et alj2002h . and increases with tem- 
perature as Mj oc (Tigm/m) 3 ^ 2 ' where \i is the mean molecular 
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weight. Using this simple scaling with temperature we obtain 

M,,(T) = M Jj0 



, ^ 3/2 / \ -3/2 



V TlGM.O 



3.8 x 1O 9 M (^|) 



Tigm\ 3/2 / n h 2 



0.15 



-1/2 



(3) 



n b h 2 V 3/5 /i + z\ 3 / 2 



0.0223 



10 



for the ionized IGM. 

The actual filter mass differs somewhat from this instanta- 
neous Jeans mass since the mass scale on which baryons succeed 
in collapsing out of the IGM along with the dark matter must be 
determined, even in linear theory, by integrating the differential 
equation for perturbation growth over time for the evolv ing IGM 
dShapiro et al1ll994l ; iGnedin & Hui|[l998l ; iGnedid 1200(1 . In real- 
ity, determining the minimum mass necessary for a halo collaps- 
ing inside an ionized and heated region to acquire its fair share of 
baryons which subsequently cool further to form stars is even more 
complicated. It depends on the detailed, non-linear, gas dynamics 
of the process and on radiative cooling. There is no single mass 
above which a collapsing halo retains all its gas, and below which 
the gas does not collapse with the dark matter. Instead, simulations 
show that the cooled g as fraction in halos decreases gradually with 
decreasing halo mass tefstathioulll992l; iThoul & Weinberdfl996l ; 
iNavarro & Steinmetzlll997l ; lDiikstraet alj|2004) . The typical halo 
sizes at which thi s transition occurs as deri ved by these different 
studies also vary. Thoul & Weinber 5 dl996h found that photoion- 
ization suppresses star formation in halos with circular velocities 
below ~ 30kms _1 , and decreases the cooled gas mass fraction in 
larger halos, with circular velocities up to ~ 50kms _1 . The halo 
circular velocity is related to its virial temperature as follows 

^ = 3.072 xl0 4 K (^)(^) 2 , (4) 

(e.g. llliev & Shapircl200lh . where ^ — 0.59 is the mean molecular 
weight for ionized gas, and the virial temperature is related to the 
halo mass by 

2/3 / 2 \ 1 / 3 



Tvir = 4.17x10 



^10 9 M© J \ 0.15 J V 10 

INavarro & Steinm etz ( 1997) found that the cooled gas fraction is 
affected by photoionization even in larger galaxies, with circu- 
lar velocities up to ~ 100 — 200 kms -1 . On the other hand, 



Dijkstraet al. (2004) recently showed, using the same method as 
Thoul & Weinberg! dl996l) , that at high redshifts the suppression 



is not as effective, and somewhat smaller galaxies can still retain 
some cooled gas. For simplicity, we assume that star formation is 
suppressed in halos with masses below 1O 9 M0 and not suppressed 
in larger halos, in rough agreement with the linear Jeans mass esti- 
mate for 10 4 K gas and the above dynamical studies. 

2.2 Source efficiencies and Pop. Ill to Pop. II transition 

We model the sources by assuming a constant mass-to-light ratio. 
Each halo found in the simulation volume at a given time, which is 
not suppressed by Jeans mass filtering is a source. For a source with 
halo mass M and lifetime t s we assign ionizing photon emissivity 
according to 



N-. 



* K fJU ( M V 7M!V (l±i\ 

VO.59/ ^10 9 M Q / ^0.15 J V 10 7 



(5), 



where the proportionality coefficient / 7 reflects the ionizing pho- 
ton production efficiency of the stars per stellar atom, Ni, the star 
formation efficiency, /* , and the escape fraction, / csc : 

A = f*fescNi. (7) 

(e.g. lHaiman & Holder! I2003T) . All these quan tities are st i ll quite 
uncertain, especially at high redshift, see e.g. Illiev"et al. (2005a) 
for discussion. Recent theoretical studies have indicated that the 
first, metal-free sta r s (Pop. Ill) mig ht have been quite massive 
(Bromm et al. 2002; Abel et al. 2000). Massive stars are more ef- 
ficient producers of ionizing photo ns, emitting up to Nj ~ 10 5 
ionizing photons per stellar atom dBromm et alj |200"II ; ISchaerej 
120021 ; FVenkatesa n & Truranl 120031) . Integrating over a top-heavy 
IMF for Pop. Ill stars leads to estimates of Ni ~ 25, 000 - 90, 000 
dSchaeredl2002h . The Salpeter IMF fo r Pop. II stars gives N t = 
3, 000 - 10, 000 dLeithereretai]|l999l) . The values of /* and / csc 
are even less certain, ranging from ~ 0.01 to ~ 1 for each of these 
quantities. There are also indications that the photon escape frac- 
tion is mass-dependent and significantly higher for small galax- 
ies at high redshift than for large galaxies o bserved at later times 
( Kitayama et al. |2004| ; lAlvarez e t al. 2006 a). Thus, many of the 
currently viable reionization scenarios involve an early population 
of small sources with high ionizing efficiency, which eventually 
evolve into the population of less efficient emitters we see at later 
times. In this work we adopt / 7 = 2000 (corresponding to e.g. 
Ni = 50, 000, /» = 0.2 and / csc = 0.2, i.e. top-heavy IMF and 
relatively efficient star formation and photon escape) for modelling 
the high-efficiency emitters and / 7 = 250 (corresponding to e.g. 
Ni = 25, 000, /* = 0.1 and / csc = 0.1, or Ni = 6, 000, /» = 0.2 
and /csc = 0.21, i.e. either moderately top-heavy IMF and moder- 
ate efficiencies, or Salpeter IMF and relatively high, but not unrea- 
sonable, efficiencies). 

The detailed mechanisms of this possible transition from high 
to low efficiency emitters, and even the question of whether it 
actually occurred are still unclear. We therefore have simulated 
cases both with and without such a transition. In the latter case, all 
sources, both large and small, have the same ionizing photon pro- 
duction efficiency. The possible physical mechanisms which result 
in decreasing the ionizing efficiency are quite varied and complex, 
including e.g. production, expulsion and mixing of metals, which 



(6) 



modifies the stellar IMF, and thus Ni, and the increase of the mean 
halo mass over time as a consequence of hierarchical structure for- 
mation, which is expected to decrease the escape fraction. Mod- 
elling and studying the detailed features and mechanisms of this 
transition is not the aim of this paper. Instead of trying to model 
these complex processes and the many related uncertainties, we 
adopt a simple model which should nonetheless capture the main 
features of the efficiency transition, as follows. For the simulations 
with varying ionizing efficiency, we assign the high value of / 7 
to the sources smaller than some characteristic mass, and the low 
efficiency to the larger halos. Since the CDM structure formation 
proceeds hierarchically, the small halos form first, and gradually 
merge up to form ever larger halos. This process of merging is pre- 
sumed to lead to a gradual increase in the mean metallicity of the 
halo gas, and decrease of the escape fractions. Physically, it is ex- 
pected that t his transition is inhomoge neous in space and extended 
in time (e.g. iFurlanetto & L oeb 2005) and proceeded faster in the 
high-density peaks, which are the first to form halos and are the first 
sites of vigorous merging. We assume the characteristic halo mass 
at which the efficiency transition occurs to be 10 9 Mq. For com- 
putational simplicity, we chose it to be the same one as the Jeans 
filtering mass since this procedure yields only two types of sources 
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(small, efficient and suppressible and large, less efficient and un- 
suppressible), rather than multiple source types. This assumption 
does not affect any of our main qualitative conclusions, though the 
exact value of the characteristic mass may affect some of the de- 
tailed quantitative estimates. If this transition mass were somewhat 
lower, e.g. 5 x 1O 8 M , rather than 10 9 Af Q , then the largest of the 
suppressed halos (between 5 x 10 8 M© and 10 9 Mq) would have 
Pop. II efficiency, while the smaller ones would have Pop. Ill effi- 
ciency. This would simply yield a case somewhere in-between our 
cases of high and low Pop. Ill efficiencies and our main conclusions 
would hold, although some quantitative numbers would change. 

2.3 Sub-grid gas clumping 

We also study the effect of gas clumping at very small scales, be- 
low the resolution of our current simulations. This clumping would 
increase the recombination rate. The sub-grid clumping coefficient, 

Csub-grid = (n 2 )/ (n) 2 , we use is given by 

n I \ IT /ice -0.114z+0. 001328 z 2 

Csub-grid {z) = 27.466e . (8) 

for WMAP1 cosmology (a good fit for 8 < z < 40) and 

r< ( \ oc om rr -0.1822z+0. 003505 z 2 , s 

C S ub-grid(z) = 26.2917e . (9) 

for WMAP3 cosmology (a good fit for 6 < z < 30), in which case 
the structures form later. These fits to the small-scale clu mping fac- 
tor are a more precise version of the one we presented in llliev et al.l 
(2005a). To derive it we used a PMFAST simulation with the same 
computational mesh, 3248 3 , and number of particles, 1624 3 , but 
a much smaller computational volume, (3.5 h _1 Mpc) 3 , and thus 
much higher resolution. These parameters correspond to particle 
mass of 1O 3 M , minimum resolved halo mass of 10 Mq, and a 
spatial resolution of ~ 1 kpc comoving. This box size was cho- 
sen so as to resolve the scales most relevant to the clumping - on 
smaller scales the gas would be Jeans smoothed, while on larger 
scales the density fluctuations are already present in our computa- 
tional density fields and should not be included again. 

The expressions in equations <[8j and (|9) exclude the matter 
residing inside collapsed halos since these contribute to the recom- 
bination rate differently from the unshielded IGM. The minihalos 
are self-shielded, which results in their lower contribution to the to- 
tal number of recombin ations than one wou l d infer from a sim ple 
gas clumping argument dShapiro et alj]2004L llliev et alj|2005bl) . In 
principle, the additional consumption of ionizing photons by mini- 
halos can also b e included as a s ub-grid prescription, as we have 
done elsewhere dCiardi et alj|2006h . which results in a further de- 
lay of the final overlap. As discussed by llliev et alJ j2005d) . how- 
ever, the biased clustering of minihaloes around the larger mass 
source halos tends to make the minihalo photon consumption cor- 
rection degenerate with the efficiency for source-halos to release 
ionizing photons, so we will assume here that this efficiency pa- 
rameter takes approximate account of this effect. The larger halos, 
on the other hand, are assumed to be ionizing sources, and their re- 
combinations are implicitly included in the photon production ef- 
ficiency, through their escape fraction. We can neglect the gas 
density fluctuations associated with suppressed sources. This gas 
is presumed to have been prevented by pressure forces from col- 
lapsing out of the IGM into these suppressed halos, so it should not 
contribute to the clumping factor as if it were inside the halos. Pres- 
sure forces may also have affected the clumping of the diffuse IGM 
inside the H II regions. Our estimate of the clumping factor, there- 
fore, serves to bracket the effect of IGM clumping. A fully self- 



Table 2. Simulation parameters and global reionization history results for 
runs with WMAP3 cosmological parameters. Box sizes are in [ /i~ 1 Mpc]. 
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consistent treatment requires following the detailed gas dynamics 
coupled to the radiative transfer, but this would require numerical 
resolution ( £ 1 kpc resolution at ~ 100 Mpc scale) which is cur- 
rently unfeasible. 

2.4 Simulation Cases 

In this paper we present a total of 15 simulations, of which 1 1 uti- 
lize the WMAP1 background cosmology. Five of these have a sim- 
ulation box size of 100 ft -1 Mpc and were presented in Papers I and 
II. This corresponds to a particle mass of 2.5 x 10 7 M Q and mini- 
mum resolved halo mass of 2.5 x 10 9 Mq (requiring 100 particles 
or more to make sure halos are properly identified). The other six 
simulations have a smaller box size of 35 h' 1 Mpc. This allows for 
better mass and spatial resolution, corresponding to a particle mass 
of 10 6 Mq and minimum resolved halo mass of 10 s Mq. This 
last mass roughly corresponds to the (generally redshift-dependent) 
minimum halo mass required for its gas to be able to cool efficiently 
by hydrogen line cool ing (a gas virial temperature of ~ 10 4 K) (e.g. 
llliev & Shapiroll200lh . These simulations and their basic parame- 
ters and features are summarized in Table [fl Note that two of our 
simulations, f2000_406 and f2000_250S_406 have the higher grid 
resolution of 406 3 for the radiative transfer calculation but are oth- 
erwise identical to simulations f2000 and f2000_250S, respectively, 
and thus serve to study possible resolution effects. 

Additionally, we investigate the effects of varying the back- 
ground cosmological parameters by doing four additional simula- 
tions for which we adopt the WMAP3 background cosmology but 
otherwise make the exact same assumptions about the reionization 
sources as the corresponding WMAP1 simulations. These simula- 
tions also impose periodic boundary conditions on the ray-tracing, 
rather than the transmissive ones adopted for the simulations in Ta- 
ble Q] The periodic boundary conditions in the radiative transfer 

2 The integrated electron-scattering optical depths for our large-b ox simu- 
lations were estimated incorrectly in Uiev et al. (2006b) and Melle ma et alJ 
l2006bl) . the correct, slightly lower, values are listed in Table[T] 
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Table 1. Simulation parameters and global reionization history results for simulations with WMAP1 cosmology parameters. Box sizes are in [ h _1 Mpc]. 
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simulations are implemented by (logically) positioning each ioniz- 
ing source at the center of the grid using the periodicity of the den- 
sity field, before calculating its contribution to the global ionization 
rates. This results in each source having a region of influence of the 
same size as our box. Since our boxes are large, they are optically- 
thick along most lines-of-sight throughout most of the evolution 
and thus the majority of the radiation is absorbed within the box. 
Any radiation that still leaves the computational volume is collected 
and put back in as a diffuse background, which boosts the effective 
local photoionization rates. The simulations using WMAP3 back- 
ground cosmology are summarized in Table[2] 

The radiative transfer simulations presented in this work were 
run on a variety of computers at CITA and The University of 
Texas, from single- and dual-processor Opteron 64-bit worksta- 
tions (at effective 3.6 GHz), to quad-processor Itanium-2 server (at 
1.3 GHz), to a 32-processor DEC-Alpha machine (at 733 MHz). 
The run times varied from 500 up to ~ 15, 000 Opteron-equivalent 
processor-hours for the simulations with 203 3 mesh, and about 8 
times longer than that for the simulations with 406 3 mesh. 

The maximum number of sources which we ray-trace per 
time-step for the large-box [(100 h _1 Mpc) 3 volume] simulations 
is up to ~ 3 x 10 5 for WMAP1 cosmological parameters, and up to 
~ 10 5 for WMAP3 cosmological parameters. The corresponding 
numbers for the smaller-box, [(35 h^ 1 Mpc) 3 volume] simulations 
are somewhat lower, at ~ 1 x 10 5 for WMAP1 and ~lx 10 4 for 
WMAP3. The net cumulative source-halo episodes explicitly ray- 
traced over the course of each simulation is typically a few million 
for the large-box simulations, and a few hundred thousand for the 
smaller-box simulations. 



3 A SIMPLE ANALYTICAL TOY MODEL FOR 
SELF-REGULATED REIONIZATION 

The effect of the suppression of small-mass sources when their ha- 
los form inside an H II region can be illustrated by a simple, an- 
alytical toy model. This will anticipate the phenomenon of self- 
regulation of reionization and the effects of varying the photon re- 



lease efficiencies of large and small halos and the IGM clumping 
factor on the evolution of the mean ionized fraction of the uni- 
verse. The universe is assumed to be comprised of H II regions 
which fill a fraction x v of the total volume and H I regions in 
the remaining volume. The global average rate of change of the 
ionized volume fraction with time, dx v /dt, is determined by the 
rate of emission of ionizing photons and the rate of recombinations 
(all per atom in the universe). The volume-averaged rates of emis- 
sion of ionizing photons per collapsed atom can be different for 
low-mass and high-mass source haloes, respectively, if the photon 
release efficiencies, / 7i i and / 7 .2 per atom are different. In addi- 
tion, the low-mass source haloes are subject to Jeans-mass filtering, 
which means that the only atoms which should be counted as col- 
lapsed onto the low-mass haloes are those associated with haloes 
in the neutral regions, which occupy a volume fraction of the uni- 
verse equal to 1 — x v . la the simplest approximation, in which the 
low-mass haloes are assumed to be distributed uniformly or ran- 
domly in space, the volume-averaged ionization rate contributed 
by the low-mass haloes is (1 — x v )fi(t), if fi(t) is the emission 
rate per atom if we neglect suppression. Recombinations occur only 
in the ionized regions, so the volume-averaged recombination rate 
per atom is the one for fully-ionized gas, fsit), multiplied by x v , 
where fs(t) = t^c — C(z)tih&b- This yields a simple differen- 
tial equation for the mean reionization history of the IGM, 

^ = Mt)(i -x v ) + f 2 {t) - f 3 (t) Xv , (io) 

where fa (t) is the emission rate of the high-mass haloes per atom 
in the universe. Equation U0\ is conveniently rewritten in a non- 
dimensional form as 

dx 
dy 

where 

c _ /i(t)+/a(t) 



-x + S, 



(11) 



h{t)+h{ty 



and 



dy = [fi(t) + Mt)]dt. 



(12) 



(13) 
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The formal solution of equation (TTJ is given by 

x v (y) = x v (0)e- y + f V dy'e-^-^S(y'), 
Jo 

where 

■t 

\h{t') + h{t')]dt'. 



(14) 



(15) 



This solution is familiar since equation dl It is identical to the static 
equation of radiative transfer in which x v is intensity, y is opti- 
cal depth and S is the source function. The contribution of the 
/i (t) term to the integral in equation dl4t is the cumulative num- 
ber of ionizing photons released by all the low-mass sources over 
time if suppression by Jeans-mass filtering is neglected, £0,1 = 
J fi(t')dt', while the integral over /3(f) is just the cumulative to- 
tal number of recombinations per atom in a fully-ionized universe, 
a large number. As long as we start the time integration at some fi- 
nite cosmic time after the Big Bang (e.g. the recombination epoch), 
this recombination integral remains finite. In that case, y ^> 1 and 
x v (0) — 0, and we can replace x v (t) by S[t), yielding 



x v (t) = 



hit) + hit) 



(16) 



h{t)+Hty 

According to equation d!6| >, reionization ends at overlap epoch 
t ov such that x v = 1, when f2(t ov ) = h(t nv ). while equation dToT> 
tells us that (dx v /dt) ov — 0. 

The solution without suppression is given by first setting 
fi(t) = in the equations above and then replacing f2(t) by 
hit) + hit): 

hit) + hit) 



[xvit)]no supp. 



hit) 



(17) 



In this case, reionization ends (i.e. x v — 1) when fi (t) + ji (t) = 
hit), and the derivative dx v /dt is again zero. 

A comparison of the solutions above for the cases with and 
without suppression shows several things. Suppression delays the 
completion of reionization, since the emission of ionizing photons 
by the low-mass sources adds to that of the high-mass haloes to 
reionize the universe and balance the recombinations at an earlier 
epoch when suppression is neglected. With suppression, the over- 
lap epoch is determined entirely by the balance between the emis- 
sion contribution of the high-mass haloes and recombinations, in- 
dependent of the contribution from the low-mass haloes. This il- 
lustrates just what is meant by "self-regulation", since the overlap 
epoch in the presence of suppression is independent of the effi- 
ciency for photon release by the low-mass haloes. Apparently, if 
they are more efficient at releasing photons, they are also more 
efficient at suppressing the formation of other low-mass sources, 
so their net effect is the same. In either case, overlap is achieved 
by the high-mass sources, which are free of suppression. Prior to 
overlap, however, the volume ionized fraction x v is higher at ev- 
ery epoch when suppression is neglected. If we were to neglect the 
small-mass sources altogether, we would find that reionization is 
also delayed, of course. But these results indicate that if we add the 
contribution of small-mass sources but account for their suppres- 
sion, the resulting overlap will be the same as if we neglected the 
small-mass sources altogether. 

Reionization is extended by the presence of low-mass haloes, 
however, relative to the case with no low-mass haloes. According to 
equations l l 1 6t and (T7), x v is higher at all times when /1 7^ than 
when /1 = 0, because ih + h)/ih+h) > h/h, since / 3 > f 2 
at all times prior to the instant of overlap (i.e. since x v < 1 before 



overlap, when x v — 1). This means that the effect of adding low- 
mass haloes would be to increase the integrated electron-scattering 
optical depth r cs without changing the overlap redshift. 

Similar simple analytical models, based on equation llOt . 
with some variations, has been explored in several recent works, 
albeit none derive d the analytical solution presented in equa- 
tions (fl6t and <fl7l jHaiman & Holdej2003l:IWvithe & Loebkool 
lOnken & Miralda-Escudel |2004 iFurlanetto & Loebll2005h . These 
studies largely reached analogous conclusions to ours, namely that 
the Jeans-mass suppression extends reionization, but only rarely, 
if ever, does this lead to a non-monotonic evolution of the ionized 
fraction, or to a double reionization. 

In reality, the reionization history is more complicated than 
our toy model suggests. As we shall see, the source haloes and 
their H II regions are spatially-clustered and thus biased relative to 
the matter distribution, so our assumptions above of uniformly dis- 
tributed low-mass haloes and uniform gas density are not correct. 
We will revisit this issue in § [4] and the Appendix and show how 
we can improve the toy model with hindsight from our detailed 
simulation results. None of the above previously-published semi- 
analytical reionization models above include the effects of bias. 



4 SIMULATION RESULTS 

4.1 Mean reionization history milestones 

We start by examining the mean global reionization histories de- 
rived from our simulations. These can be characterized by several 
basic parameters. The first of these parameters is the epoch of over- 
lap, z ov , which we define as the time when the mass-weighted ion- 
ized fraction of the gas first surpasses 99%. This also quantifies 
the overall duration of reionization, since its start is determined by 
when the first resolved halos form in our simulations, and thus is 
fixed by the structure formation alone. The second global parame- 
ter is the total electron scattering optical depth, r cs , integrated from 
the beginning of reionization to the present. The third and final pa- 
rameter is the redshift when 50% of the gas mass is ionized for the 
first time. This is of direct interest for observations since this half- 
ionization point is a good indicator of the epoch when the fluctua- 
tions of the redsh ifted 21 -cm emission from neutral hydrogen reach 
their maximum dMellema et alj2006bl) . 

For all the large-box simulations of comoving size 
100 h" 1 Mpc on a side the first resolved halo forms, and thus reion- 
ization starts, at z ~ 20 (16) in WMAP1 (WMAP3) cosmology. 
For the smaller-box simulations of size 35 /i -1 Mpc this occurs 
earlier, at 2 ~ 30 for WMAP1 cosmology and at z ~ 22 for 
WMAP3 cosmology. This earlier start for smaller simulation vol- 
ume is due to the hierarchical nature of CDM structure formation, 
whereby the smaller halos form earlier. Our N-body simulations 
have the same number of particles, thus the smaller boxes have cor- 
respondingly higher mass resolution, and hence the first halos form 
earlier. Among our simulations, the earliest epoch of overlap is 
z ov = 13.5, achieved in the (somewhat artificial) case f2000_250, 
i.e. when small-mass sources are never suppressed by Jeans mass 
filtering, while at the same time they are highly efficient at produc- 
ing ionizing photons. The combination of these two properties re- 
sults in a very large photon emissivity, and as a consequence, a very 
fast reionization and too early overlap. In the more realistic cases, 
in which the small-mass source suppression is included, the overlap 
occurs significantly later, between redshift z — 10.4 (f2000_250S) 
and 8.4 (f250C_250S). For the large-volume simulations, the red- 
shifts of overlap are similar, ranging from z — 11.3 (f2000) to 
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Figure 1. Mean reionization histories: (a)(left) 100 hT x Mpc box runs: (bottom panel) redshift evolution of the mass-weighted ionized fraction, x m for f2000 
(solid, blue), f2000C (dotted, red), f250 (short-dashed, green), and f250C (long-dashed, magenta), (top panel) Corresponding ratios of mass-weighted and 
volume-weighted ionized fractions, which are equal to the mean density of the ionized regions in units of the mean density of the universe, and (b)(right) 
35 ft -1 Mpc box runs: 0000.250 (solid, blue), f2000.250S (dotted, red), f2000C.250S (short-dashed, green), f250.250S (long-dashed, magenta), and 
f250C_250S (dot-short dashed, cyan). Insets: the same reionization histories, but in linear scale. 



z = 8.2 (f250C). The reason for these similar end-of-reionization 
times despite the lack of small-mass sources in the large-volume 
simulations is that most of these small sources become suppressed 
by the time of overlap, and thus, in either case, the reionization 
is brought to an end by the same large-mass sources. However, 
at early times there are significant differences between the results 
for the two box sizes. As we noted above, the small-mass sources 
start forming significantly earlier, thus models which resolve these 
sources have higher ionized fraction at early times. This results in 
higher integrated optical depths when small-mass sources are re- 
solved. For the WMAP1 cases, these range from r cs = 0.173 
(f2000_250) when the small-mass sources are high-efficiency but 
not suppressed, to r cs = 0.148 (f2000_250S) when suppression 
is taken into account, to r cs = 0.111 (f250C_250S) when low- 
mass sources have the same efficiency as the high-mass ones and 
sub-grid IGM clumping is taken into account, all of which are 
roughly within l-a from the WMAP1 estimate, r cs = 0.17 ± 0.04. 
In the 100 h^ 1 Mpc box simulations, the small-mass sources are 
not present and thus reionization starts later. This results in some- 
what lower values of the integrated optical depth, ranging from 
t cs = 0.130 (f2000), which is nevertheless still in agreement with 
the WMAP1 measurement, to r cs = 0.098 (f250C), which is a bit 
low. 

Increasing the resolution to 406 3 grid has only a modest effect 
on the reionization histories. In these cases the underlying density 
field is resolved better, which results in an increased recombination 
rate, and a slightly more extended reionization history. 



4.2 WMAP1 vs. WMAP3 and ionizing source efficiencies 

Within the WMAP3 background cosmology, structure forma- 
tion occurs later and all the evolution shifts to correspondingly 
lower redshifts. For the same efficiencies, overlap is now reached 
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Figure 2. Mean reionization histories: 35 h 1 Mpc box runs with WMAP3 
background cosmology: f2000.250S (dotted, red), f2000C.250S (short- 
dashed, green), and f'250_250S (long-dashed, magenta). Inset: the same 
reionization histories, but in linear scale. 



at redshifts between z ov = 7.9 (f2000_250S) and z ov ~ 7 
(f2000C_250S), in rough agreement with the data from high- 
redshift galaxies and SDSS QS O's, which indicate a tail-end of 
reionization at z QV ~ 6 - 7 (e.g. lFan et al.l2002l ; lwhite et al.l2003l; 
Malh otra & Rho ads 2004). The resulting integrated electron scat- 
tering optical depths are correspondingly lower, at r cs = 0.082 — 
0.103, in complete agreement with the WMAP3 constraints. These 
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results confirm that the delay in structure formation due to the lower 
values of erg and Qo, higher value of h and the slight tilt of the pri- 
mordial power spectrum found by WMAP3 results yield a decrease 
of the integrated optical depth from the WMAP1 value of ~ 0.17 
to the WM AP3 value of ~ 0.0 9, as was previously predicted an- 
alytically dAlvarez et a"il20 06b). This analytical estimate predicted 
that the change of the cosmological parameters from the WMAP1 
to WMAP3 values would result in a delay of reionization by a fac- 
tor of 1.4 in (1 + 2) and a corresponding lower optical depth by 
a factor of 1.4 3 / 2 ~ 1.7. These predictions are nicely confirmed 
by our current simulations, which yield a decrease in the redshift 
of 50% ionization by a factor of ~ 1 .3- 1.4, and in the overlap red- 
shift by a factor of ~1.3, while the r cs values decrease by factors of 
~ 1.4 ~ 1.3 3 ' 2 . The correction factors are slightly smaller than the 
ones derived analytically. These differences are due mostly to the 
use of periodic boundary conditions in our WMAP3 simulations 
(which ensures that no photons are lost through the simulation box 
boundary), vs. our use of transmissive boundary conditions in the 
earlier simulations. This yields a slightly faster evolution at the late 
stages of reionization, and correspondin gly earlie r overlap. 

There has been a recent claim (Popa i2006h that the delay in 
structure formation due to the different background cosmological 
parameters derived by WMAP3 is not in fact sufficient to account 
for the reduction of the derived optical depth, once the feedback 
and radiative transfer effects are accounted for. This work claimed, 
instead, that the WMAP3 cosmology requires a more top-heavy 
IMF in order to match the new value of the integrated optical depth. 
Our simulations show conclusively that this is not the case. 

Our results show that, while a very top-heavy IMF is not re- 
quired, ionizing sources nonetheless must have been fairly efficient 
photon producers. This means that, compared to the present-day 
observed galaxies, the high-redshift galaxies must have had either 
a moderately top-heavy IMF, a higher star formation efficiency, a 
higher photon escape frac t ion, o r some combination of these. 

Recently. IZahn et alj 12006) performed a simulation similar to 
the ones we presented in Papers I and II. This work used an ap- 
proach similar to ours in coupling a radiative transfer scheme to the 
results of an N-body simulation of the density field. They adopted 
WMAP1 cosmology parameters, but utilized a smaller simulation 
volume and a different radiative transfer method. The sources re- 
solved in these simulations have minimum mass 2 x 10 9 M Q , very 
similar to the resolution of our large-box simulations. In terms of 
ionizing source efficiencies, they assumed that all the sources only 
produce a cumulative one photon per every atom in the universe by 
redshift z ~ 6.5. This corresponds to a much lower source emis- 
sivity than the ones we have assumed. Such a model predicts quite 
late reionization and a very low integrated electron scattering op- 
tical depth of r es = 0.06. The final overlap is achieved by z ~ 6 
only if the number of recombinations per atom integrated over time 
is negligible. As we noted above, adopting the WMAP3 cosmology 
for these simulations would inevitably push the redshift of overlap 
down to ~ 4 and the optical depth down to 0.035 (and even fur- 
ther down if recombinations are included), in clear disagreement 
with the observations. This again demonstrates that relatively high 
photon production efficiencies are required for the high redshift 
sources, similar to the ones we adopt here. 



4.3 Globally-averaged reionization histories 

The full reionization histories (mass-weighted ionized fractions vs. 
redshift) from all of our simulations are shown in Figs. [T] and [2] 
The character of reionization is similar in all cases, regardless of the 
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Figure 3. Effect of Jeans-mass filtering on the collapsed mass fraction of 
halos in our simulations. Shown are the cumulative collapsed mass fractions 
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collapsed fraction at that time) vs. M for all source halos (solid, black), and 
the unsuppressed source halos for cases f2000_250S (short-dashed, red) and 
f'250_250S (dotted, blue) at three redshifts, as labelled. 



simulation volume or the background cosmology. In all cases reion- 
ization is clearly inside-out, with the high-density regions being 
ionized on average earlier than the low-density regions, as shown 
by the ratio of the mass-weighted to the volume-weighted ionized 
fraction, x m /x v , which reflects the mean density of the ionized 
bubbles comp ared to the mean density of the universe, as reported 
previously in Jlliev et al.ll2 006b). These ratios are plotted in the top 
panels of Figs. [TJ and [2] In all cases and at all times the ionized 
regions are overdense, more significantly so at early times. This 
inside-out character of the reionization process is more pronounced 
for the small-box simulations which have higher spatial resolution 
and thus follow the underlying density field more closely. 

The evolution is always monotonic, and in no case we do see 
a double reionization or even a modest temporary decrease of the 
ionized fraction. The evolving source efficiencies (from a high one 
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Figure 4. Collapsed fractions in low-mass halos (10 8 ^ M/Mq ^ 10 9 ; 
dotted) and high-mass halos (M/Mq > 10 9 ; solid) for the simulation with 
(35 h _1 Mpc) 3 volume and WMAP1 cosmology. 

early-on to a low one later), small-scale gas clumping (from low 
values early to higher values later), and the Jeans-mass filtering of 
sources all only slow the process down and extend it in time, but 
never manage to reverse it even temporarily. The high-resolution, 
small-box simulations find somewhat more extended reionization 
than the large-box simulations, due to the earlier onset of reioniza- 
tion in these cases, combined with the suppression of small-mass 
sources. 



4.4 Self-regulation 

One interesting and important consequence of the low-mass source 
suppression is that the time of overlap becomes fairly insensitive 
to the efficiency of the small-mass sources, since these are almost 
completely suppressed by the reionization end. Furthermore, in 
cases when the small-mass sources are more efficient photon emit- 
ters, their H II regions expand quickly, thereby suppressing fur- 
ther source formation in a larger volume than less efficient sources 
do. This is illustrated in Figure [3] where we show the cumulative 
source halo collapsed fraction with and without suppression, nor- 
malized to the corresponding totals, at three different redshifts - 
early, middle and late in the evolution. During the early evolution, 
the high-efficiency of the small sources helps suppress a signifi- 
cantly larger fraction of them, by a factor of a few or more, com- 
pared to the low-efficiency case. At late times, the level of suppres- 
sion becomes very large and essentially independent of the small- 
source efficiency. By then most of the volume is ionized, and due to 
the strong bias in their spatial clustering, concentrating them in the 
ionized regions, the effect on the small-mass sources is even more 
dramatic than simply the ionized volume would suggest. This sup- 
pression notably slows down the further expansion of the ionized 
bubbles, thus resulting in a very efficient self-regulation of reion- 
ization and a relative insensitivity of the global evolution to the 



assumed small-mass source efficiency. For example, if we com- 
pare the cases f2000.250S and f250.250S (with WMAP1 param- 
eters), we find that the ratio of their mass ionized fractions starts 
at ~ 7 — 8 early-on (z > 20), decreasing to ~ 3 at z ~ 15 and 
to ~ 1.5 or less for z < 13. The assumed small-mass source effi- 
ciency is 8 times larger in the former than in the latter case, while 
the larger, unsuppressible sources have the same efficiency in the 
two cases. This shows that during the very early evolution the ratio 
of the ionized fractions matches the ratio of the number of photons 
produced, but then the ionized fractions ratio quickly decreases as 
the evolution progresses and the low-mass source suppression and 
the corresponding self-regulation start to dominate the reionization 
process. With the WMAP3 background cosmology we observe the 
same generic behaviour, except that it is shifted to later times, as 
we discussed above. 

In Figure|4]we show the total collapsed fractions in low-mass 
and high-mass halos for the 35 h~ Mpc computational box with 
WMAP1 parameters. At early times the collapsed fraction in small- 
mass halos rises faster and without suppression they dominate the 
large-mass halos at all times. The collapsed fraction in high-mass 
halos rises exponentially below z — 15. The collapsed fractions in 
the two mass ranges become comparable at z ~ 8. When suppres- 
sion is included, large-mass halos begin to dominate the reioniza- 
tion progress much earlier, after z ~ 14 — 15, depending on the 
small-mass source efficiency, as shown in Figure[3]above. 

In Figure [5] we show the evolution of the cumulative number 
of ionizing photons, £, emitted by all sources in our 35/i _1 Mpc 
simulation volume, for both the WMAP1 and the WMAP3 cases. 
The suppression of sources due to radiative feedback clearly has 
a quite dramatic effect on the cumulative £, and can reduce it by 
up to 2 orders of magnitude in some cases compared to the case 
when no sources get suppressed. The effect is stronger when the 
small-mass sources are more efficient at producing ionizing pho- 
tons, again demonstrating the self-regulated nature of reionization 
- more efficient sources suppress more, partially compensating for 
the higher efficiency. The emissivity per unit time and per baryon 
in the universe is given by d£ / dt. This derivative is strictly positive 
at all times since d£/dt = (d£/dz)(dz/dt), dz/dt < and based 
on Fig.[5]also d^/dz < 0), so the emissivity per baryon rises with 
time in all models we consider here . 

In Figure|6]we show the cumulative number of recombinations 
per atom in our computational volume. The number of recombina- 
tions is initially small since a very small fraction of the volume 
is ionized. This is despite the fact that these first H II regions are 
highly-overdense, and thus their recombination times are short. At 
later times, the recombinations become more important, consuming 
0.5-1 additional photon per atom when sub-grid gas clumping is 
ignored (which underestimates the recombinations) and up to 5-10 
additional photons per atom when sub-grid clumping is included. 
This demonstrates that recombinations play an important role dur- 
ing reionization and should not be ignored in any simulations or 
analytical models. 

Some of these trends for the basic self-regulation of reioniza- 
tion were roughly anticipated by our analytical toy model in § [3] 
That model neglects several effects which our detailed simulations 
show are quite important, however. These include the inhomoge- 
neous gas density inside the ionized regions and the clustered for- 
mation of source haloes in time and space, which biases their distri- 
bution and that of the H II regions relative to the underlying matter 
density field. In order to account for some of these neglected effects 
in an approximate way, to gain further insight and use the toy model 
quantitatively, we have incorporated the results of our simulations 
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Figure 5. The cumulative number of ionizing photons per total atom, f, emitted by all sources in the 35 ft -1 Mpc box simulations with WMAP1 (left) and 
WMAP3 cosmology (right). Same notation as in Figure[T](b) and Figure|2] respectively. For reference we also show the total £ for f'250_250S and f25()C_250S 
simulations if no sources were suppressed (dot-long dashed, black). 
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Figure 6. Cumulative number of recombinations vs. redshift for: (a)(left) 100 h 1 Mpc box runs, f2000 (solid, blue), f2000C (dotted, red), f250 (short- 
dashed, green), and f250C (long-dashed, magenta), and (b)(right) 35 ft" 1 Mpc box runs, f2000.250 (solid, blue), f2000.250S (dotted, red), f2000C.250S 
(short-dashed, green), f250.250S (long-dashed, magenta), and f250C.250S (dot-short dashed, cyan). 



to improve the simple analytical model described in § [5] (see Ap- 
pendix for details). In this improved model the volume- weighted 
ionized fraction x v evolves according to: 

^ = 6i(t)/i(t)(l -xv) + fr(t) - f 3 (t)x v , (18) 

where bi(t) = (1 — a;™)/(l — x v ) is the correction factor 
which accounts for the fact that low-mass halos are concen- 
trated toward the ionized regions rather than being distributed uni- 
formly, or randomly in space. The photon emission rate fi (t) = 



ft, small (/coil, small/ A£) is the mean photon emissivity per atom 
coming from the small-mass sources based on our modelling of 
sources in the radiative transfer simulations, where At is the time 
between two time slices (~ 20 Myr), over which we hold the un- 
derlying gas density field and halo population fixed while we inte- 
grate the radiative transfer and the non-equilibrium ionization bal- 
ance rate equations. In the toy model in § [3] the source suppres- 
sion by Jeans-mass filtering is assumed to be limited to the ionized 
volume, with source halos randomly, or uniformly, distributed in 
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Figure 8. Spatial slices of the ionized and neutral gas density from simulations (a)(left) f2000_406 (100 h Mpc comoving) at redshift z = 12.9, and 
(b)(right) from f'2000_250 (35 h^ 1 Mpc comoving) at z = 16.2. Both are at similar ionization fraction by mass of x m ~ 60%. Shown are the density field 
(green) overlayed with the ionized fraction (red/orange). The cosmic web of structures is visible in both the neutral and ionized regions. 




Figure 9. Spatial slices of the ionized and neutral gas density from simulations (a)(left) f2000_250S_406 (35 h^ 1 Mpc comoving) at redshift z = 14.1, and 
(b)(right) from f2000_250S_406 (35 h^ 1 Mpc comoving) with WMAP3 parameters at z = 10.1. As in Figure[8] these images are also at ionization fraction 
by mass of x m ~ 60%. Shown are the density field (green) overlayed with the ionized fraction (red/orange). 

space, thus the unsuppressed fraction is proportional to the neutral 
volume fraction, 1 — x v (i.e. n = 1 in eq. flOl ). We also con- 
sider a second case, where we empirically model the effects of halo 
spatial clustering bias by setting n — 0.1 in equation dlOt (see Ap- 
pendix). The second term, /2(f) = /-y,large(/coil/At) is the same 
as /1 (t) but for the larger, unsuppressible sources. Finally, the last 
term, /3(f) = l/t ICC — C(z)nHttH is the rate of recombinations 
per atom (assuming fully-ionized gas within the H II regions at the 
mean density). Recombinations only occur in the ionized volume, 
and so the volume-averaged rate is just proportional to x v . In re- 



ality, the H II regions are overdense on average, which we ignore 
here, and thus our model somewhat underestimates the effect of 
recombinations. 

Solutions for the evolution of the mean ionized fraction based 
on our toy models are plotted in Figure [7] for the small-box sim- 
ulation cases in Table [T] These curves show that when low-mass 
sources and their suppression are included, but bias is neglected 
(n = 1), the end of reionization is nearly the same, regardless of 
the efficiency adopted for photon release per low-mass halo atom 
(see curves for cases f2000_250S and f250_250S). When bias is in- 
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Figure 7. Reionization histories based on a simple analytical model which 
incorporates sources with different efficiencies, source suppression, and 
evolving gas clumping. Notation is the same as in Fig.[T](right panel). Thin 
lines correspond to the cases with empirical "bias" included (see Appendix 
for details). 

corporated (n — 0.1) there is a very modest trend toward an earlier 
end of reionization when the low-mass source efficiency is higher. 
When bias is neglected, the reionization evolution tends to be too 
fast compared to the simulations; bias serves to slow reionization 
down because a larger fraction of the low-mass sources are sup- 
pressed at any given time than the neutral volume fraction 1 — x v , 
since these halos are concentrated in the H II regions. The effect of 
gas clumping is to delay the end of reionization, but since clump- 
ing grows with time and the contribution from recombinations in 
the ionized regions is small at early times, soon after the sources 
turn on, the evolution at early times is similar with and without 
clumping. 

4.5 Effect of small-mass sources on the reionization 
morphology 

The morphology of the H II regions during reionization shows 
many common features, but also some interesting variations 
(Figs.[8]and[9]l. In general, H II regions persist and grow over time, 
as new sources form to add their ionizing luminosity to the total 
emission rate inside each ionized volume. In all cases reionization 
occurs inside-out, i.e. the denser regions around the high-density 
peaks are ionized first. Locally, at any given scale the ionized bub- 
bles percolate at very different times. This local overlap leads to 
quite large (of size 10 Mpc comoving, or larger) ionized regions 
fairly early, well before the whole universe is ionized. This explains 
why previous, small-box simulations failed to see the gradual build- 
up of the ionized volume of the universe over time, prior to the fi- 
nal overlap epoch, and mistakenly found that reionization is a very 
rapid transition. At all times prior to overlap, there are also many 
isolated sources or small groups of sources which create a number 
of small H II regions. The low-density regions remain largely neu- 
tral during most of the evolution, since they are devoid of sources. 
The case when only the large (and unsuppressible) sources are 



present (case f2000_406; Figure[8] left) has been discussed in some 
detail in Paper I. When smaller sources are also present, the result- 
ing geometry depends on the radiative feedback upon them. When 
the low-mass sources are not suppressed, the basic morphology is 
similar to that when only larger sources were present (Figure [8}, 
with some additional small-scale features due to the now-resolved 
small-mass sources. However, when the small-source suppression 
is included, the morphology changes, even with the same source 
efficiencies (Figure|9]l. There is more small-scale structure present, 
although the large, locally-percolated regions are of similar sizes 
to the ionized regions observed in the large-box simulations, indi- 
cating that by this time reionization is already dominated by the 
large-mass sources. In the presence of the small-mass, suppress- 
ible sources, there are also small, "relic" H II regions occasionally, 
i.e. regions which were ionized earlier but are now recombining be- 
cause some (or all) of the sources which made them originally were 
suppressed (these relic H II regions are seen as darker spots, mostly 
around the edges of the large ionized bubbles). Generally, these 
relic H II regions are short-lived, however, being quickly overrun 
by the I-fronts from neighbouring ionized bubbles, and thus do not 
have strong effects on the evolution. 

A direct comparison of the morphology of reionization in the 
simulations presented here with tho se by other meth ods, such as 
reported by jCiardi et alj|2003bl) or jZahn et alj [2006) is difficult, 
since some or all of the properties which determined their outcome, 
including box size, resolution, source efficiencies and cosmologi- 
cal parameters, are generally different. Of these, moreover, only 
the simulations presented here were able to resolve the low-mass 
source halos subject to suppression and took this suppression into 
account. We shall present a more detailed discussion of the char- 
acteristic scales and topology of reionization based upon our sim- 
ulation results in a forthcoming paper, along with further compar- 
isons. The purpose of the current paper, however, is to demonstrate 
the differences made by improving the resolution to take account 
of leans-mass filtering, as described above. 



5 CONCLUSIONS 

We have studied the effects of small-mass sources on the progress 
and duration of reionization based on a large suite of detailed radia- 
tive transfer simulations. We found that these small-mass sources 
play an important role during reionization. In their presence reion- 
ization starts much earlier, by Az ~ 10 in WMAP1 cosmology (at 
z ~ 30 rather than at z ~ 20), and by Az ~ 6 in WMAP3 cosmol- 
ogy (at z ~ 22 rather than at z ~ 16). They also supply most or all 
of the ionizing emissivity at early times, when the larger galaxies 
are still exceedingly rare. 

However, the same low-mass sources are also a subject to sup- 
pression in the ionized regions due to Jeans-mass filtering. This 
low-mass source suppression decreases the global ionizing photon 
emissivity by factors of up to 100, depending on the photon pro- 
duction efficiency of the small-mass sources. As a consequence, 
the Jeans-mass filtering of small-mass sources effects the global 
reionization history profoundly, delaying overlap by Az ~ 3, and 
decreasing the integrated electron-scattering optical depth r cs by 
~ 0.025 compared to the case without suppression. The evolv- 
ing gas clumping at small (sub-grid) scales extends the reioniza- 
tion significantly, by Az ~ 1 — 1.5, but has only a modest effect 
on the corresponding electron scattering optical depth, decreasing 
it by ~ 0.013 - 0.014. 

Furthermore, the reionization process is strongly self- 
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regulated - the more efficient the small-mass sources are, the larger 
the fraction of them that are suppressed. This results in a corre- 
sponding decrease in the number of ionizing photons emitted, par- 
tially canceling the effects of the high efficiency. However, the later 
stages of reionization are completely dominated by the larger-mass 
sources which do not get suppressed. Their numbers rise exponen- 
tially as they become more common, while, in contrast, the low- 
mass halos are strongly suppressed at late times. As a result, the 
overlap epoch is not very sensitive to the properties of the low- 
mass sources, but instead is largely determined by the numbers and 
efficiencies of the high-mass sources. As a consequence, for a given 
overlap redshift there is a significant degeneracy between the small- 
mass source's efficiency and the r cs - by changing the efficiency 
of the small-mass sources we can easily increase or decrease the 
global integrated optical depth without changing the overlap red- 
shift. 

In short, the end of reionization is dictated by the higher-mass 
haloes which are not subject to Jeans-mass filtering, but the addi- 
tion of small-mass, suppressible haloes serves to extend the reion- 
ization epoch to earlier times, thereby boosting r cs . This naturally 
explains why WMAP CMB polarization measurements indicate a 
large r cs and early beginning for reionization, z > 1 1 for WMAP3, 
while the absorption spectra for quasars at z 6 indicate that the 
overlap epoch was at z < 7. Previous simulations failed to see this 
because they failed to include the small-mass halo sources with 
Jeans-mass filtering in a simulation volume big enough to sample 
the global reionization history fairly. 

The combined effects of low-mass source suppression in the 
ionized regions, decreasing ionizing source efficiency and increas- 
ing gas clumping over time can cause reionization to "linger" at a 
given ionization fraction for some time, before eventually reaching 
overlap. Ho wever, none of our mod els s how the do uble reionization 
predicted bv lWvithe & Loebl J2003I) and lCerj ( l2(5o3l) . What prevents 
double reionization from occurring is the large spread in the reion- 
ization histories from one region to another, combined with the 
strong source clustering and the self-regulation. Double reioniza- 
tion scenarios require that the highly-efficient Pop. Ill sources emit 
enough photons to reionize the universe at high redshift, but then 
die or become less efficient, causing the universe to partially re- 
combine and be finally reionized later by the rise of Pop. II sources. 
However, these highly-efficient sources form at very different times 
in different places, and thus the transition from high-efficiency to 
low-efficiency sources is spread over time. These early sources also 
get suppressed quite efficiently, thus never get to large enough total 
emissivity at any given time, so as to reionize the whole universe. 
These self-regulation mechanisms should hold under quite generic 
circumstances, beyond our s pecific realizations. Some of these 
were previously discussed bvlFurlanetto & Loebl j2005h within a 
simplified analytical framework, which did not include some im- 
portant effects like the strongly increased Jeans-mass filtering due 
to source clustering. Nevertheless, they reached conclusions similar 
to ours, namely that double reionization is not physically plausible. 

Reionization proceeds in an inside-out fashion, whereby the 
high-density peaks b ecome ionized first and voids last, thereby con- 
firming the results o f llliev et al.l < f2006bl) . The mean overdensity of 
the ionized regions is always greater than one, and up to 4 at early 
times. This trend is even more pronounced in the smaller-box, high- 
resolution simulations, where the density fluctuations are resolved 
better. In this case also the shapes of the H II regions are more non- 
spherical during their early evolution. Since at early times both the 
ionized density field and the source distribution are strongly bi- 
ased relative to the underlying density field, it is crucial to perform 



detailed numerical simulations in order to obtain the correct ge- 
ometry, size and spatial distributions of the ionized bubbles. Most 
current semi-analytical models of reionization ignore the effects of 
halo bias on the Jeans-mass filtering, which leads to a large un- 
derestimate of the low-mass sou rce suppression. The only excep- 
tion is the very recent work of Kra mer et al.l d2006l) which pre- 
sented a model with an approximate, spherically-averaged source 
bias model. While such models cannot give the full H II region 
geometry, it would be interesting to compare their results to full 
simulations in order to see to what extend their statistical results 
are representing the bias effects reliably. 

The overlap by redshift z ~ 6 — 7 indicated by current obser- 
vations is easily achieved by stellar sources with either a Salpeter 
or a slightly top-heavy IMF. There is, therefore, no need for addi- 
tional, more "exotic" sources of radiation (e.g. decaying particles) 
in order to satisfy both the electron scattering optical depth con- 
straint from WMAP and the end-of-reionization constraints com- 
ing from the SDSS QSO's and high-redshift Ly-a surveys, as they 
currently stand. 

In WMAP3 cosmology the formation of cosmic structures, 
and thus the epoch of reionization, are delayed compared to 
WMAP1 cosmology, due to a combination of the lower normal- 
ization of the power spectrum of density fluctuations, tilt, lower 
matter density and slightly higher Hubble constant. This delay is 
such that it almost exactly compensates for the lower value of r es 
found by WMAP3, given ionizing sources with the same efficiency. 
This was predicted analytically bv lAlvarez et al.l y006bj) and con- 
firmed by our current simulations. 

Our results show that our large-box simulations which include 
only the large-mass, unsuppressible sources, nevertheless predict 
correctly the overlap epoch and the large-scale features of reion- 
ization. This feature is a consequence of the strong self-regulation 
of reionization, due to which by the time H II regions grow large 
(roughly when x v £ 0.1) the vast majority of the low-mass sources 
are already suppressed and the reionization process is driven by the 
large-mass sources. However, any simulations which do not include 
the low-mass sources would underestimate the ionized fraction, 
particularly at early times, and the small-scale features of reion- 
ization, both of which are influenced by the low-mass sources. The 
total integrated electron-scattering optical depth is thus also under- 
estimated. Semi-analytical reionization models which do not in- 
clude the self-regulation properly would similarly underestimate 
the mean ionized fraction history and the corresponding optical 
depth. On the other hand, the large-scale features of reionization are 
not very sensitive to the low-mass sources and their efficiencies, but 
depend significantly on the large-scale features of the density field, 
the clustering of the ionizing sources and level of dumpiness of 
the gas, all of which are most reliably modelled through large-scale 
reionization simulations. 
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Figure Al. The evolution of the fraction of low-mass ionizing sources 
which are Jeans-mass suppressed based on our simulations (solid lines, as 
labelled), a simple model, where the suppressed fraction is simply the mean 
volume-weighted ionized fraction (dotted lines), and an empirical model 
where this fraction is x®' 1 , to account for source bias (dashed lines). 
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APPENDIX A: A SIMPLE TOY MODEL FOR 
REIONIZATION WITH SELF-REGULATION AND HALO 
BIAS 

In order to gain some further insight into the self-regulation of 
reionization and its dependence on the assumed parameters, we 
construct the following simple toy model. The global average rate 
of change of the ionized volume fraction with time, dx v /dt, is 
determined by the rate of emission of ionizing photons and the 
rate of recombinations (per atom in the universe). The ionizing 
source model we use in our radiative transfer simulations assumes 
that their emissivity is proportional to the collapsed gas fraction 
in halos, / co n, with a given efficiency, / 7 (see § 12.21 . In this, we 
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distinguish the sources whose formation is suppressed inside the 
H II regions by Jeans-mass filtering from the high-mass sources, 
whose formation is not suppressed. In this framework the rates 
of emission of ionizing photons from the small-mass and large- 
mass sources are given by /i = (/ 7 ) sm aii/coii, S maii/At and 
f'2 = (/ 7 ) large /coil, large/At, respectively, where At is the time 
over which these photons are emitted. The former rate should also 
be modulated by the ionized fraction, since small sources are sup- 
pressed in ionized regions. In the simplest approximation, one can 
assume that the suppressed fraction of small-mass sources is sim- 
ply proportional to the ionized volume fraction, x v , i.e. the emit- 
ting (unsuppressed) fraction is given by /i(l — x v ). Such a model 
would be exact if the small-mass sources were randomly distributed 
in space. In practice, however, the halos that host the sources are 
generally clustered together, thus the suppression of small-mass 
sources should be stronger than a random distribution would give, 
as these sources cluster around the high-density peaks, which are 
ionized earlier in the inside-out reionization yielded by simulations. 
Therefore, in an effort to crudely model the effects of this source 
bias, we also consider a second case, where we modulate the sup- 
pression by a factor of &i(t) = (1 — — x v ), where we set 
n = 0.1 (n = 1 corresponds to the random distribution above). 
This power of 0.1 was determined empirically, by roughly match- 
ing the actual Jeans-mass suppressed fraction from our simulations 
to 1 — x". In Figure lATl we plot both models against the actual data 
from the simulations (for WMAP1 parameters). We note that this 
simple model is not intended to match the simulation results per- 
fectly, but rather to show how the mean reionization histories vary 
under different assumptions. Thus, even though it does not follow 
the Jeans-suppressed source fraction exactly, it it is appropriate for 
our purposes here, as it reflects the main evolutionary trends seen 
in our simulations. 

Finally, the rate of recombinations per unit time per atom is 
given simply by the inverse of the recombination time, = t^ ec = 
C(z)riHCiB, where C(z) is the (evolving) clumping factor, as is 
the Case B recombination rate of hydrogen, and we assume full ion- 
ization in the H II regions. Combining these factors, we can write 
the equation for the evolution of the volume-weighted ionized frac- 
tion x v as: 

^ = /i(t)(l-xJ) + / a (t)-/ 8 (t)i„ (Al) 

Based on our simulation data (in WMAP1 cosmology), the 
functions fi and fa as function of redshift are well-fit by the fol- 
lowing expressions: 

fi =0.335 ex P (0.227z - 0.02546z 2 ) Myr" 1 (A2) 

h =4.6312 ( ^ 7 2 ' 5 q SC ) ex P (-0.107* - 0.02463z 2 ) Myr _1 ,(A3) 
while fz is given by 

* = l^W^- CA4) 

6.655 x 10 a 

Finally, in order to use these expressions directly, we need to 
change the independent variable in equation I A 1 1 from time to red- 
shift, using 

dy dydz dy . . , 3 ,i /2 .... 

dt = TzTt =-^o(l + z)[fio(l + *) +V A ]' , (A5) 

as appropriate for flat cosmology. The results given by this model 
are plotted in Figure|7]and discussed in the main text. 



